##########################################################################################
# Title: Figure 3a — Population Density within 5 km of a Business District
# Section: 2 (Figures)
# Purpose: Bar chart by city for 2000/2010/2020; sorted by 2020 (desc), exclude Mumbai
# IO:
#   IN : Processed_Data/Intermediate/Section_2/Processed_BD_Analysis/
#        ├─ ALL_CITIES_100K_09_2000.rds
#        ├─ ALL_CITIES_100K_09_2010.rds
#        └─ ALL_CITIES_100K_09_2020.rds
#   OUT: Graphs/Figure_3a.png
##########################################################################################

rm(list = ls()); options(stringsAsFactors = FALSE)

suppressPackageStartupMessages({
  library(data.table)
  library(ggplot2)
  library(grid)   # unit()
})

# ---- Project root (align with run_all.R) -------------------------------------
if (!requireNamespace("here", quietly = TRUE)) install.packages("here", quiet = TRUE)
ROOT      <- normalizePath(here::here(), winslash = "/")
PDIR      <- file.path(ROOT, "Processed_Data", "Intermediate", "Section_2")
IN_DIR    <- file.path(PDIR, "Processed_BD_Analysis")
GRAPH_DIR <- file.path(ROOT, "Graphs")
dir.create(GRAPH_DIR, recursive = TRUE, showWarnings = FALSE)

# ---- Inputs ------------------------------------------------------------------
f2000 <- file.path(IN_DIR, "ALL_CITIES_100K_09_2000.rds")
f2010 <- file.path(IN_DIR, "ALL_CITIES_100K_09_2010.rds")
f2020 <- file.path(IN_DIR, "ALL_CITIES_100K_09_2020.rds")
stopifnot(file.exists(f2000), file.exists(f2010), file.exists(f2020))

out_file <- file.path(GRAPH_DIR, "Figure_3a.png")

# ---- Helper: extract population_d_5 and convert to thousands/km^2 ------------
pull_pd5 <- function(path, year_label) {
  dt <- readRDS(path) |> as.data.table()
  dt[bounds == 5L & !is.na(city_name) & !is.na(population_d_5),
     .(city_name, value_thous = population_d_5 / 10, year = year_label)]
}

d2000 <- pull_pd5(f2000, "2000")
d2010 <- pull_pd5(f2010, "2010")
d2020 <- pull_pd5(f2020, "2020")
df    <- rbindlist(list(d2000, d2010, d2020), use.names = TRUE)

# ---- Filter + reorder --------------------------------------------------------
df <- df[city_name != "Mumbai" & !is.na(value_thous)]
# sort by 2020 descending, then reverse for top-to-bottom ordering after coord_flip
order_2020 <- df[year == "2020"][order(-value_thous), city_name]
df[, city_name := factor(city_name, levels = rev(order_2020))]  # <-- key change
df[, year := factor(year, levels = c("2000","2010","2020"))]

# ---- Colors ------------------------------------------------------------------
cols <- c("2000"="#B0C4DE", "2010"="#4682B4", "2020"="#0B3C6F")

# ---- Axis + reference line ---------------------------------------------------
xmax  <- max(df$value_thous, na.rm = TRUE)
upper <- ceiling(xmax / 5) * 5
ref_vals <- c(10)

# ---- Plot --------------------------------------------------------------------
p <- ggplot(df, aes(x = city_name, y = value_thous)) +
  geom_col(aes(fill = year), width = 0.70, position = position_dodge(width = 0.80)) +
  coord_flip(clip = "off") +  # prevents right-edge truncation after flip
  scale_fill_manual(
    values = cols, name = NULL,
    guide = guide_legend(direction = "horizontal", nrow = 1,
                         label.position = "bottom",
                         keywidth = 3.0, keyheight = 0.35, label.hjust = 0.5)
  ) +
  labs(
    y = "Population Density within 5 km of a Business District (Thousand ppl/km²)",
    x = NULL
  ) +
  scale_y_continuous(
    limits = c(0, upper),
    breaks = seq(0, upper, by = 5),
    expand = expansion(mult = c(0, 0.12))  # extra headroom on the right (post-flip)
  ) +
  geom_hline(yintercept = ref_vals, linetype = "dotted", color = "grey60") +
  theme_classic(base_size = 10) +
  theme(
    legend.position      = "bottom",
    legend.justification = "center",
    legend.box           = "vertical",
    legend.box.spacing   = unit(0.25, "lines"),
    legend.text          = element_text(margin = margin(r = 18)),
    plot.margin          = margin(t = 18, r = 48, b = 26, l = 10),  # more room on right
    axis.title.x         = element_text(margin = margin(t = 6)),    # nudge title from edge
    panel.background     = element_rect(fill = "white", color = NA),
    plot.background      = element_rect(fill = "white", color = NA)
  )

ggsave(out_file, p, width = 6, height = 4, dpi = 300, bg = "white")
message("✅ Wrote: ", normalizePath(out_file, winslash = "/"))
